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■i \ Abstract - We present a new method for the numerical treatment of 

c : 

O . second order phase transitions using the level spacing distribution function 

O 

P{s). We show that the quantities introduced originally for the shape analysis 
of eigenvectors can be properly applied for the description of the eigenvalues 
as well. The position of the metal-insulator transition (MIT) of the three 
dimensional Anderson model and the critical exponent are evaluated. The 
shape analysis of P(s) obtained numerically shows that near the MIT P{s) is 
clearly different from both the Brody distribution and from Izrailev's formula, 
and the best description is of the form P{s) = ci sexp(— C2 s^"*"^), with (3 ~ 0.2. 
This is in good agreement with recent analytical results. 



I. INTRODUCTION 

Recently a novel method has been introduced for the location of the crit- 
ical point and the determination of the critical exponent in the three dimen- 
sional (3D) Anderson model exhibiting a metal-insulator transition (MIT). It 
has been demonstrated by Shklovskii et al} and by Hofstetter and Schreiber^''^ 
that random matrix theory (RMT)^ may serve as a tool for a surprisingly ac- 
curate calculation. This time the necessary information is derived using the 
spectrum rather than the wave functions of the system. It is also expected 
that this new method is easily applicable for other types of second order phase 
transitions. 

The novel approach is based on the study of the statistical properties 
of the eigenvalues of the Hamiltonian on both sides of the MIT. Based on an 
analogy between the kicked rotator and the Anderson model^ the MIT can be 
considered as a transition from the chaotic regime to the non-chaotic one, or in 
other words using the terminology of the RMT, from the Gaussian orthogonal 
ensemble (GOE) to the Poisson ensemble (PE) of random matrices. The level 
spacing distributions for both ensembles are known: for the GOE it is very well 
described by the Wigner surmise 

PGOE(s) = ^sexp(-^s2), (1) 

that shows linear level repulsion for low s. For the PE one has 

PpE(s)=exp(-s), (2) 

i.e. in this case the energy levels of localized states may be arbitrarily close. 

The model under consideration is described by the usual tight-binding 
Hamiltonian 

l-L = ^ei\i><i\+^V\i><j\, (3) 



where i labels the sites of a simple cubic M x M x M lattice. In the second sum 
only the nearest neighbor interactions are considered, and for sake of simplicity 
we chose V^ = 1 as the unit of the energy scale. The potentials Ei are the site 
energies taken from a uniform distribution —W/2 < Si < W/2. Therefore the 
disorder W will be the critical parameter. 

One expects that for small disorder the spectrum of the Hamiltonian 
should be described by the level statistics of the GOE where due to hybridiza- 
tion level repulsion occurs and states become delocalized, while for large enough 
disorder the eigenvalues will tend to be uncorrelated random numbers and the 
corresponding eigenstates will be localized. Therefore as disorder increases the 
MIT is accompanied with a transition from Pgoe{s) to Ppe(s) with some un- 
known spacing distribution Pce{s) at the MIT. (The index CE stands for the 
critical ensemble occuring, as demonstrated below, at the MIT.) Pce{s) may at 
the same time show characteristics of both the GOE and the PE as suggested 
by Shklowskii et al.^ In infinite systems this transition is discontinuous;^ how- 
ever, simulations in finite systems show a continuous variation of the level 
spacing distribution. In fact there is a scaling property^"^ of these P{s) as 
M changes for any fixed value of W. The sign of this scaling is clearly seen 
as a fixed ensemble in the -P(s), namely the CE obtained for different values 
of M. Moreover, there appears a fixed point sq ~ 2 in the P{s) curves for 
different disorders W. Therefore one may divide the interval [0, oo) into [0, 2] 
and [2, oo). The first part has been studied in Refs. 2, 3 and the latter (which 
is equivalent due to the normalization of P{s)) in Ref. 1. This time we will use 
all the numerically obtained P{s) functions over a wide interval s G [0, 5]. 

The transition between the GOE and the PE can be approximated by 
several interpolation formulas. One of them is due to Brody^ 

PB{s) = c^s^eM-C2s'+^), (4) 
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where ci and ci are determined according to the conditions of normahzation 

P(s) ds - 1 (5) 

and that the mean spacing is unity 

/•OO 

/ sP{s) ds = 1. (6) 

Jo 

Any spacing distribution P{s) should satisfy Eqs. (5) and (6). Therefore we 
have C2 = [T{{/3 + 2)/(/3 + 1))]^+^ and ci = (1 + /3) C2. Another interpolation 
formula was given by Izrailev^ 

Pi{s) = As^{l + Bl3sfe-^''-^', (7) 

where 

and the constants A and B are to be calculated numerically according to con- 
ditions (5) and (6). Both of these interpolations give back the limiting cases: 
for (3 = 1 the GOE distribution and for /? = the PE one. 

Concerning the P{s) close to the MIT, in a recent publication Aronov et 
al.^ have shown analytically that the distribution at the transition may well be 
described by 

PA(s) = cisexp(-C2s'+^), (9) 

where constants ci and C2 are fixed according to conditions (5) and (6) and 
for parameter (3 they obtained < /3 < 1. Furthermore fi is related to the 
correlation length exponent u by^ 



In this paper we wish to show a numerical analysis that confirms the form of 
Eq. (9) and at the same time provides a critical exponent that satisfies relation 
(10). 

II. SHAPE ANALYSIS OF THE LEVEL SPACING DISTRIBUTION 

Since we expect to see a transition from the GOE to the PE statistics as 
disorder increases we propose to study such quantities that describe the shape 
of the calculated P{s) and compare them to the known limiting cases. If one 
parameter scaling holds the plot of these quantities versus disorder obtained for 
different system sizes M should show a fixed point, yielding the approximate 
position of the critical point as well as the approximate value of the critical 
exponent. For such a calculation Shklovskii et al} used the tail, s G [2, oo) of 
P(s), while Hofstetter and Schreiber^ employed the numerical fit of Eq. (7) 
with Eq. (8) on the other part, s G [0, 2], of the P{s). The latter authors have 
also analyzed^ the integrated level statistics and the Dyson-Mehta statistics, 
as well, and shown that these quantities enable an even better finite size scaling 
then P{s). In this contribution we introduce a different approach for the char- 
acterization of the level statistics. We will use all of our numerically obtained 
P{s) functions and in contrast to previous methods we will not introduce spe- 
cial parameters other than those that are uniquely related to the shape of the 
distribution function of a set of random numbers. 

It has already been shown^'^ that it is advantageous to characterize a 
set of non-negative random numbers by certain moments of their distribution. 
This problem may arise studying e.g. noisy wave functions. The quantities 
introduced in Ref. 10 are the spatial filling factor or participation ratio which 
is calculated as 

At2 



and the structural entropy 

5,,, = ^+ln^, (12) 

Hi Hi 

where fxi and H2 are the usual first and second moments of the distribution 
p{x) of the random variables 

POO 

l^k = I x^ p{x)dx, (13a) 

Jo 

and Hs is calculated as 

/•OO 

Hs = — xln{x) p{x)dx. (13b) 

Jo 

As Eq. (6) ensures /Ui = 1 when using the level spacing distribution we will 

have simply q = ^2^ a-i^d Sstr = /Ug + ln(U2. Note that in practical calculations 

Eqs. (13a) and (13b) are approximated by finite sums. The shape analysis 

resides on the comparison of points plotted on the (g, Sstr) plane with curves 

calculated with known p[x) functions. ^^ Note that for a trivial distribution 

p{x) = d{x—xo), one obtains q = I and Sstr = 0. For any other distribution one 

will have q < 1 and Sstr > and the relations < q < 1 and < Sstr < — In g 

always hold. 

The above characteristics will be employed here as well. Eqs. (11) and 
(12) should be calculated for every value of disorder parameter W and system 
size M replacing x by s and p{x) by P{s). The calculated values (q^Sstr) 
will be compared to the continuous curves obtained using the interpolation 
formulas due to Brody (4) and to Izrailev (7) as well as with other possible 
P{s) functions. We will show that the Pce{s) is qualitatively different from 
the ones obtained for GOE and PE. 

We would like to emphasize that the calculation of q and Sgtr is a method 
which is not affected by the position of the fixed point in the P{s) at s ~ 2. 
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Therefore we are rather using aU the obtained level distribution. At the same 
time no fitting procedure is necessary. We should note that one of the quanti- 
ties, q, originally used as the participation ratio of a wave function serves as a 
measure of the skewness (or peakedness) of the distribution in our context. The 
structural entropy has no previously known meaning in our formalism. Addi- 
tionally we have to mention that these quantities contain information about 
many-level correlations.^ 

III. ANALYTICAL CALCULATIONS 

First we give the (g, Sstr) relations for the interpolating distributions as we 
wish to compare the numerical results with these phenomenological functions. 
As /3 runs from zero to unity the application of definitions (9) and (10) using 
Pb('S) from Eq. (4) yields explicitly 



.B 



?"(/3) 






2 



n'f^l. 



,B ,__^i'(3 + 3\ ^^f(3 + 2\ 1 _J(3 + 2 



and 

Here ifj^x) is the digamma function. On the other hand the q^{(3) and 5'gj^(/3) 
functions obtained using Pi{s) from Eqs. (7) and (8) can be calculated numer- 
ically with sufficient accuracy. Both of the cases are shown in Fig. 1, where 
the quantities, q and Sstr are plotted as a function of /?. 

As one can observe following the transition from PE to GOE the form 
of the P{s) changes in two ways: first the low-s behavior changes from a 
constant to linear level repulsion and at the same time the large-s tail changes 
from exp(— s) to exp(— s^). These two changes are accounted for by both of 
the interpolation functions (4) and (7). However, at the transition in our 
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physical system one might expect that Pce{s) shows characteristics of both of 
the two Umiting ensembles. Keeping this in mind we introduce further possible 
interpolations between the exponential -PpeIs) and the Wigner Pqoe{s), e.g. 
such an intermediate distribution (IMl) may look like 

PiMi{s) =cis^exp{-C2s), (16) 

where ci — 1 + (3 and C2 — ci/r{l + (3). The parameter /? runs in the [0, 1] 
interval. This distribution is that of the PE for /3 = and at /3 = 1 it has the 
low-s behavior of the GOE. Similarly another intermediate distribution (IM2) 
may look like (see also Eq. (9)) 

PiM2{s)=cisexp{-C2S^+f^), (17) 

where ci = (l+/3)[r(3/(l+/3))]V[r(2/(l+/3))]3 and cs = [r(3/(l+/3))/r(2/(l+ 
/9))]^+^. In this formula the parameter (3 also runs in the [0,1] interval. At 
(3 = 1 this distribution is just the Wigner surmise and for (3 = it coincides 
with PiMi('S) with /3 = 1, i.e. the two functions meet with 

PIM(s)=4se-2^ (18) 

which for s ^ 1 is a GOE and for s ^ 1 is a PE distribution. The quantities 
q and Sstr for these new intermediate P{s) distribution as a function of their 
parameters are 

^'""'m = 1^, (19) 



for PiMi(s) and 



^ir(/3)=ln(/3 + 2)-^(/3 + 2) (20) 



^ ^^ ' r(2/(i + /3))r(4/(i + /3))' ^ ' 



.-„).. „r(^)-,nr(^)-^.(^) (22) 

for PiM2{s). We have plotted Eqs. (19), (20) and Eqs. (21), (22) in Fig. 1. 
The combination of the two intermediate forms can be given as a two-parameter 
form 

PiM3(s)=cis''exp(-C2s"), (23) 

but for sake of simphcity we restrict ourselves to the one-parameter versions of 
either Eq. (16) or (17). 

IV. NUMERICAL CALCULATIONS AND RESULTS 

In our investigation we have used the results of the numerical simulation 
presented and described in detail in Ref. 2. We have taken the data of the 
P{s) histograms and calculated the quantities q and Sstr as a function of W 
around the critical disorder 15 < VF < 18 for different system size M ranging 
from 13 up to 21. The states were obtained at the band center (E = 0) for 
which the critical disorder is expected^^ to be around Wc = 16.5. 

First we present our results concerning the position of the critical point 
and the critical exponent. In Fig. 2 we show for M = 21 how the calculated q 
and Sstr values change with the increase of disorder interpolating between the 
PE and GOE values. For an infinite system one expects a step function-like 
behavior, here it is smeared out by the finite size of the system. In Fig. 3 
we have plotted our results for both quantities for different system sizes. The 
dashed line in both figures shows the expected position of the critical disorder. 
It is clear that a fixed point exists around Wc = 16.75 for both quantities. In 
this work we have calculated the fixed point from second order polynomial fits 
to the data and averaged over the different pairs of M and M' . The value one 
obtains in this way is Wc = 16.87 ± 0.52 for q and Wc = 16.77 ± 0.63 for Sstr- 



The critical exponent can be determined in a similar way. The approxi- 
mate value is given by 

_ InjM/M') 
'^^'^' - ln(AM/AMO' ^ ^ 



where 

dX 



A 



M 



dW 



(25) 



M,M' 



W^ M' is the approximate value for the critical disorder obtained for the pair of 
M and M' . The averaged results for the critical exponents are u = 1.27 ± 0.29 
for X = q and u = 1.30 ± 0.38 for X = Sstr- 

From these calculations we may conclude that this method does indeed 
give quantitatively correct results. The value of z/ ^ 1.34 for the critical ex- 
ponent is obtained in Ref. 3. Note that the resulting value for the critical 
disorder Wc is slightly higher than 16.5 as in recent calculations of the multi- 
fractal properties of the wave functions at the MIT.^^ But noting the accuracy 
in both cases, we point out that Wc = 16.5 is well within the error bars. 

Now we analyze the calculated Sstr values as a function of q a-sW changes 
around the critical point. Such a relation may be compared with the continuous 
ones obtained in the previous section. In Fig. 4 we display the data in the Sstr 
vs q diagram. The symbols denote the simulation and the curves the analytical 
results. The numerical data show a remarkable trend in that respect that they 
fall onto a common line independent of the system size for a wide range of 
disorder 15<VF<18. It is also clear that this trend is different from the 
SstriQ) curve observed for the Brody or the Izrailev distribution. In Fig. 5 we 
have enlarged the most important part of Fig. 4. As the numerical simulation 
leads to (g, Sstr) values close to that of the approximation IM2 we conclude that 
the empirical -P(s) function should have very similar properties as the Pim2{s) 
has. In fact choosing the calculated P{s) function for M = 21 and W = 16.75, 
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the corresponding quantities are q ~ 0.703 and Sstr ~ 0.156. These values 
can be obtained with the choice of /3 = 0.18 ± 0.02 in Eq. (17), while for 
W = 16.5 we compute q ~ 0.708 and Sstr ~ 0.153 which can be reproduced 
with P = 0.21 ± 0.02 in Eq. (17). Hence we conclude that the intermediate 
distribution Pim2{s) with a parameter (3 ~ 0.20 gives a good approximation of 
the -Pce('S) at the MIT. In Fig. 6 we show that the numerical histogram at 
W = 16.75 for M = 21 is well approximated by the distribution of the form of 
Eq. (17) with (3 = 0.2. 

We note that this distribution shows the GOE characteristics for small 
level spacing s. For large s, however, does not follow the PE statistics (2) so 
that our data do not support the expectation that the CE shows characteristics 
of both limiting ensembles. 

The visible discrepancy between the curve for Pim2 and the calculated 
data in Figs. 4 and 5 is still unknown. We have performed preliminary cal- 
culations with the two-parameter distribution Pjms given in Eq. (23). This 
function can approximate the numerical points with a better accuracy, e.g. with 
the choice 6 ^ 1.3 and a ~ 1.1. Such situation with 5 > 1, however, would 
violate a general symmetry theorem by Dyson. ^^ We have also performed a 
very accurate analysis of small-s data considering the integrated level spacing 
distribution and obtained a value 6 ~ 0.97 for the best fit at the MIT.^^ 

The result presented is at the same time capable to explain the strange 
behavior of the normalization parameter A (see e.g. Eq. (7)) observed in Ref. 
2, where at the transition Hofstetter and Schreiber reported indications of a 
discontinuous change of A as a function of W. The normalization constant 
in Eq. (17) with /3 = 0.2 is ci ~ 2.28 which is larger than for the Poisson 
and Wigner distributions. In the limit M — > oo the normalization constant is 
expected to be A = 7r/2 on the metallic side and A = 1 on the insulating side 
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while at the transition it is larger than both values. Similar arguments for the 
parameter B in Eq. (7) follow from Eqs. (7) and (8) and the above arguments 
for A, explaining the respective observations in Ref. 2. 

V. CONCLUSIONS 

We have presented a general method for the analysis of the spacing dis- 
tribution around the critical point of a second order phase transition. As an 
example we have calculated the position of the MIT and the critical exponent 
in the 3D Anderson model. Although with lower accuracy, this method does 
indeed give the correct answer. On the other hand the main result of our paper 
is that for the MIT in the 3D Anderson model we have found a possible shape 
of the spacing distribution -PceIs) as given by Eq. (17). We have also pre- 
sented an explanation for the strange behavior of the normalization constant 
observed in Ref. 2. Our results are in agreement with recent theoretical expec- 
tations derived by Aronov et al.^ which yield^ the relation (10) between (i and 
the correlation length exponent v. For our numerical value z/ ^ 1.3 in (i = 3 
the relation (10) yields a value 13 = 0.26, which is close to the value (3 ~ 0.2, 
obtained from the shape analysis above. 
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Figure Captions 

Fig. 1. Quantities q and Sstr as functions of ttie free parameter f3 for tiie 
Brody (dastied line) tlie Izrailev distribution (solid line), the first in- 
termediate (IMl) (dashed-dotted line), and the second intermediate 
(IM2) distribution (dotted hue). 

Fig. 2. Calculated q and Sstr for M = 21 as a function of disorder W. 
The limiting values for the GOE and PE are shown with horizontal 
dashed-dotted lines. The expected position of the MIT is shown by 
a vertical dashed line. 

Fig. 3. Calculated Sstr (a) and q (b) for M = 13, 15, 17, 19, 21 as a function 
of disorder W. The expected position of the MIT is shown by a 
vertical dashed line. 

Fig. 4. Calculated Sstr vs q for M = 13, 15, 17, 19, 21 in the range of 15 < 
W < 18. The solid circles are the points representing the PE, the 
GOE and the IM (Eq. (18)) cases. Solid line represents the relation 
for the Izrailev, dashed line the one for the Brody distribution. The 
dashed-dotted line reflects the intermediate distribution IMl (be- 
tween PE and IM) given in Eqs. (19), (20), and IM2 (between IM 
and GOE) given in Eqs. (21), (22). 

Fig. 5. A part of Fig. 4 enlarged. The result of the simulation atW = 16.75 
for M = 21 is plotted bold and marked with the text MIT. The solid 
line is obtained from the interpolation with the Izrailev, dashed line 
with the Brody, and dashed-dotted line with the intermediate (IM2) 
distribution. 

Fig. 6. The numerically obtained histogram of P{s) for W = 16.75 and 
M = 21 (solid line) and the distribution given in Eq. (17) with 
(3 = 0.2 (dashed line). 
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